Quantum phases of mixtures of atoms and molecules on optical lattices 
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We investigate the phase diagram of a two-species Bose-Hubbard model including a conversion 
term, by which two particles from the first species can be converted into one particle of the second 
species, and vice-versa. The model can be related to ultracold atom experiments in which a Feshbach 
resonance produces long-lived bound states viewed as diatomic molecules. The model is solved 
exactly by means of Quantum Monte Carlo simulations. We show that an "inversion of population" 
occurs, depending on the parameters, where the second species becomes more numerous than the 
first species. The model also exhibits an exotic incompressible "Super-Mott" phase where the 
particles from both species can fiow with signs of superfiuidity, but without global supercurrent. 
We present two phase diagrams, one in the (chemical potential, conversion)-plane, the other in the 
(chemical potential, detuning)-plane. 
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I. INTRODUCTION 

In the past years the Bose-Hubbard model |l| has been 
extensively investigated and a lot of interest has been 
generated thanks to ultracold atom experiments on op- 
tical lattices 01, which provide an ideal realization of 
the model. Recently, much theoretical and experimental 
work has been performed on mixtures with several species 
of particles. For instance, Bose-Fermi mixtures on lat- 
tices have been studied [3|, U, l5|, d, 0, [1] ■ Another mixture 
that is likely of interest involves atoms and molecules, 
in which conversion between the two species is possi- 
ble. Such conversion processes can describe, for instance, 
long-lived bound states of atoms (diatomic molecules) oc- 
curing in ultracold atom experiments where a Feshbach 
resonance is used to tune the scattering length of the 
atoms In those experiments, the hyperfine inter- 

action between two spin polarized atoms can flip the spin 
of one of the atoms, reducing sensitively their scattering 
length. The two atoms are virtualy bound into a "molec- 
ular" state until the hyperfine interaction flips again the 
spin of one of the atoms. 



II. THE MODEL 

With the motivation above, we propose to study a 
two-boson species model with an additional conversion 
term allowing two particles from the first species to turn 
into one particle of the second species, and vice-versa. 
We denote the first species as "atoms", and the second 
species as (diatomic) " molecules" . Atoms and molecules 
can hop onto neighboring sites, interact, and conversion 
between two atoms and a molecule can occur. Several 
atoms can reside on the same site, their interaction be- 
ing described by an on-site repulsion potential. A second 
on-site repulsion potential describes the interactions be- 
tween molecules and atoms being on the same site. This 
leads us to consider the following Hamiltonian 



with 



f^-taJ2{4aj+h.c.) -i,„^(mjm^. +/I.C.) (2) 
P = UaaY^ht (ht-l) + UaraJ2n^nr+Dj2n^ (3) 



(4) 



The T, P, and C operators correspond respectively 
to the kinetic, potential, and conversion energies. The 
a| and operators (to| and m^) are the creation and 
annihilation operators of atoms (molecules) on site i, 
and "fif — ajoj (n™ = m|mj) counts the number of 
atoms (molecules) on site i. Those operators satisfy 
the usual bosonic commutation rules [a^ , aj] = 6ij and 



|m,-,m]| = S, 



jj . In order to simplify the model and re- 



H = T + P + C, 



(1) 



duce the space of parameters, we impose a hard-core 
constraint on molecules. This is done by adding the 
condition m^m^ = mjmj = 0. For a minimal model, 
we set a maximum of two atoms per site by imposing 
OjajOj — a|a|aj = 0. The sums («, j) run over pairs of 
nearest- neighboring sites i and j. We restrict our study 
to one dimension and we choose the atomic hopping pa- 
rameter ta ^ I in order to set the energy scale, while we 
choose the molecular hopping parameter tm = 1/2, mo- 
tivated by the continuous-space behavior of the hopping 
as a function of the mass {t oc /2m), a molecule be- 
ing twice heavier than an atom. Smaller values of tm (as 
mapping of experimental systems to Bose-Hubbard mod- 
els would suggest [l,[l3|) are not expected to lead to qual- 
itatively different behavior. The parameter Uaa controls 
the interaction strength between atoms, and Uam con- 
trols the interaction between atoms and molecules. The 
conversion between atoms and molecules is controlled by 
the positive parameter g. This parameter can be related 
to the "hyperfine interaction" parameter in the Fesh- 
bach resonance picture ^, Ji|. Finally, the parameter 
D acts as a chemical potential for molecules, and al- 
lows to tune the energy difference between atomic and 
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molecular states. This parameter can be related to the 
"detuning" in the Feshbach resonance example. In the 
remainder of the paper we will not expand on the connec- 
tion to the Feshbach resonance problem, nor attempt to 
reproduce Feshbach resonance physics. We concentrate 
on taking the model given in (IT])-((4]) at face value and 
determining its phase diagram. A similar model for the 
one-dimensional continuum has been analysed in Ref , 
and for optical lattices in the mean field approximation 

It is important to note that the Hamitonian ([1]) does 
not conserve the number of atoms Na — X^i'^l'^i' ^"-"^ 
the number of molecules N„i = J2i "^l^i j because of 
the conversion term However we consider that a 

molecule is made of two particles, so the total number of 
particles N in the system is conserved: 

N^Na + 2Nra. (5) 



III. QUANTUM MONTE CARLO 
SIMULATIONS: THE WORLD LINE 
ALGORITHM 

In order to make the model suitable for simulations, 
we perform a mapping of the Hamiltonian describing two 
species of bosons on a ID lattice ([TJ onto a Hamiltonian 
describing single species of bosons evolving on a ladder 
(Fig. [T|). In the ID space, the two species live together. 
They can hop onto neighboring sites, and the interac- 
tion between the two species is described by an on-site 
potential Uam- The conversion between the two species 
occurs on a single site. In the ladder space, the atoms 
(molecules) live on the top (bottom) side of the ladder. 
The interaction between the two species is described by a 
potential Uam acting between vertical neighboring sites. 
Two atoms living on the same atomic site can be de- 
stroyed at the same time, with the creation of a molecule 
on the corresponding molecular site (and vice- versa). 

• Atom • Molecule o Single species o Empty site 

Uaa Uam }^ 



I O • O )|^^ O • • O ?|; g Mixed space 



W O O O • O • • O Atomic space 

0»00*0000«''^ Molecular space 



FIG. 1: (Color online) In order to make the model suitable 
for simulations, a mapping is performed between the model 
of atoms and molecules living on a ID lattice, and a model of 
single species where the particles reside on a ladder. 

Quantum Monte Carlo simulations are performed for 
the ladder model by making use of the World Line algo- 
rithm [l^, . It is essential to emphasize that this al- 
gorithm works in the canonical ensemble, meaning here 



that the total number of particles N = Na + 2Nm is 
conserved. Indeed, simulations using a gran d canonical 
algorithm (Stochastic Series Expansion) [l^ turned out 
to be difficult to handle, because it is numerically very 
hard to control the number of particles of each species 
using two chemical potentials, the number of particles of 
each species depending on both chemical potentials. 

Defining the continuous product of evolution operators 
in imaginary time. 



M 



0-/3 



one starts by writing the partition function as the trace 
of the evolution operator e"'^^ 



dr 



^ = E (^1 n iv') 



(7) 



using the occupation number representation for the states 
[ip). Then we use the so-called "checkerboard decompo- 
sition" for the Hamiltonian, H = He + Ho, with 



He — E 



i odd 



where Hi = Ti + jPi + jCi and Ti,Pi, Ci are defined by 
Tt = -ta{a\a^^^+h.c) -tm{m\m^^^-^h.c) (9) 

P^ = Uaa [n^ {fl^ - l) + fl^^, {fl^^, - l)] (10) 

+ Uam{n1nT+n1+,n^r+i) 

+ D{nT+nT+i) 

(7, = -g{m\a^ai+a\a\m^) (11) 

- 5("^I+la^+la,+l +aJ+iaJ+iTOi+i). 

We attract here the attention of the reader to EqfTT| 
in which we have added a minus sign to the conversion 
term. The energy of the model is independent of the sign 
in (fTT|) . so (lll|) and ([4]) are equivalent. This can be seen 
by realizing that flipping the sign of the conversion term 
just results in a redefinition of the phase of the molecu- 
lar creation and annihilation operators, m]' = —m] and 



mj = —mi. We work with a minus sign in pip in order to 

ensure that all matrix elements ((^|e~'^'^|i/') are positive. 
Those positive matrix elements normalized by Z define 
the probability of transition from the state to the 
state |(/)), which is required for a Monte Carlo sampling. 

It is important to note that He and Ho are written each 
as a sum of operators Hi that commute (but He and Ho 
do not commute). Using the Trotter-Suzuki formula at 
second order, 



-dri^H^ + H^ 



drH^ 



kdrHo 



OidT^),{U) 
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and using properties of the trace we get 



dT 



(13) 



The error due to the Trotter-Suzuki decomposition 
vanishes because of the continuous product making dr 
going to zero (in the case of a discrete product the Trot- 
ter error becomes ©(dr^) instead of 0(^dT^), due to the 
accumulation of errors in the product). Introducing com- 
plete sets of states / = X]i/.(r) IV''(''')) (V'(''')l between each 
pair of exponentials leads to 



mr^T/2)) (14) 



^ = E n Wr-HiT)|e-'^-^= 

x(^(T+dr/2)|e-'^"^°|V(r)), 

where the sum runs over all sets of states ipir) for all 
values of r in [0,0]. Finally, each operator e~'^'^^= and 
^-drHo jg product of independent four-site operators 
g-drffi f2 sites i and i -f 1 in the atomic space and 2 sites 
in the molecular space). With the hard-core constraint 
on molecules and a maximum of two atoms per site, the 
size of the Hilbert space of the four-site problem is 36. 
Thus each matrix element in p4p can be computed by 
evaluating numerically 36 x 36 matrices. As a result, 
the quantum problem has been mapped onto a classical 
problem with an extra imaginary time dimension, and the 
algorithm consists in generating configurations of states 
iP{t) using standard classical Monte Carlo techniques. 
For more details, see references [T6l|. 



Indeed, this allows to keep an even total number of par- 
ticles, preventing an extra single particle to be out of the 
atoms/molecules conversion process. 

Another quantity of interest for the characterization of 
a phase is the superfluid density. An easy way to access 
this quantity is to make use of Pollock and Ceperley's for- 
mula flT'l that relates the superfluid density to the fluc- 
tuations of the winding number W, ps = L(W^') /2tf3, 
where t is the hopping of the considered species, /3 the 
inverse temperature, and L the number of lattice sites. 
Usually, this winding number W is perfectly well-defined 
for systems with n species of particles. For a given con- 
figuration, it is defined by the number of times that the 
world lines cross the boundaries of the system from the 
left to the right, minus the number of times they cross 
the boundaries from the right to the left (Fig. [2^). But 
in our case, the atomic and molecular windings, Wa and 
Wm, are ill-defined because the world lines associated to 
each of the species may be discontinous if conversions be- 
tween atoms and molecules occur (Fig.[5)D). It is then no 
longer possible to determine whether a particle is flow- 
ing to the right or to the left as a function of imagi- 
nary time. However we can define atomic and molecular 
pseudo- windings, W* and W^, by the number of right 
jumps minus the number of left jumps, normalized by 
the number of sites L. Non-zero values of such pseudo- 
windings are signatures of superfluidity of the particles. 
When no conversion between atoms and molecules oc- 
curs, the deflnition of pseudo- winding coincides with that 
of true winding. In addition, the correlated winding is 
well-defined for the mixture of particles. 



2W* 



(18) 



IV. QUANTITIES OF INTEREST 



In addition to the atomic and molecular densities, 



Pa^Na/L prn 

we also define the total density 



Ptot 



Na + 2N„ 



(15) 



(16) 



because the composite atomic and molecular world lines 
are continuous (if one considers that a molecular world 
line represents two atomic world lines). This corre- 
lated winding is relevant for the superfluid density of the 
mixture because it corresponds to the winding of par- 
ticles, without looking at their individual nature (atom 
or molecule). It is also interesting to consider the anti- 
correlated winding, 



War..t = W: - 2W1 



(19) 



by analogy with ([5]), where L is the number of sites in 
the lattice. 

In order to identify insulating phases, it is useful to 
look at the behavior of the total density ptot as a function 
of the chemical potential p{N). It is common to define 
the chemical potential in the canonical ensemble at zero 
temperature by the energy cost to add one particle to 
the system, p{N) = E{N + 1) — E{N). However, for our 
present model, it is better to define it by the energy cost 
to add successively 2 particles to the system divided by 
2, 



piN) 



E{N + 2) - EiN) 



(17) 



which allows to determine if atoms and molecules are 
flowing in opposite directions or not. The definitions of 
correlated winding (|18p and anti-correlated winding ([TJ 
are similar to those used in Bose-Fermi mixtures [6, _8J 



V. NUMERICAL RESULTS 

A. The one-site problem 

It is useful to start the investigation of the model by 
considering first the one-site problem with a total num- 
ber of particles N = 2 {ptot = 2). Figure [3] shows the 



4 



It 



It 



W=l 



W =? — Atom 

a 

W =? — Molecule 




12 3 4 
(a) 



12 3 4 



(b) 



FIG. 2: (Color online) Example of world lines for a four-site 
lattice with periodic boundary conditions, (a) For a system 
without conversion between the different species, the world 
lines are continuous and the winding number is well-defined, 
(b) The conversion between atoms and molecules leads to dis- 
continuities in the world lines, and no true winding can be 
defined for each of the species. However it is well-defined for 
the mixture atoms/molecules, because the composite world 
lines are continuous (see text for details) . 
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FIG. 3; (Color online) The one-site problem with 2 particles. 
The densities of atoms and molecules are plotted as functions 
of the conversion parameter g for differents values of the de- 
tuning D, for Uaa = 4. 



B. The lattice problem 



atomic and molecular densities as functions of the con- 
version parameter g and different values of the detuning 
D for Uaa — 4 (the value of [/am does not play any role 
since there is only 2 atoms or 1 molecule) . For D = and 
small g, the 2 particles are mainly bound in the molecular 
state, because the creation of the molecule has a vanish- 
ing energy cost while having 2 atoms costs 2Uaa = 8. As 
g increases, it becomes energetically favorable to make 
conversions atoms/molecule, so the atomic density starts 
to grow, reducing the molecular density. When g is large, 
the system maximizes the conversion process. Thus the 
system is in the molecular state with 1 molecule half of 
the time, and in the atomic state with 2 atoms the rest of 
the time. As a result, the atomic and molecular densities 
converge to pa = Ptotf^ = 1 and pm = Ptot/4: = 1/2. For 
D = 6 the same behavior holds, but the molecular density 
decreases faster to the large g limit because the energy 
associated to the molecular state is higher and closer to 
that of the atomic state. For D = 10 we have the inverse 
behavior, the molecular density increases with g and the 
atomic density decreases, because it is now cheaper en- 
ergetically to have 2 atoms rather than 1 molecule. The 
transition point between those two cases is D ~ 8 ~ 2Uaa 
for which the atomic state has exactly the same energy 
as the molecular state. Those states have the same prob- 
ability, and varying g just changes the rate of conversion 
between them. Thus the expectation values of the atomic 
and molecular densities do not depend on the value of 
and remain equal to the values that optimize the conver- 
sion process: Pa = 1 and pm = 1/2. 



We now turn to the full problem with L lattice sites. 
We have performed simulations for L — 20, 40, 80, 160 
and determined by extrapolation to L = cx) that fi- 
nite size effects associated to the choice of working with 
L — 20 lead to errors smaller than our statistical error 
bars, these latter being smaller than the size of the sym- 
bols displayed in the figures of this paper (unless other- 
wise stated). In the same manner, we have determined 
that using j3 — L allows to get the physics relevant to 
the ground state {(3 = oo), for the measured quantities. 
As for the one-site problem, we start by looking at the 
atomic and molecular densities as functions of 17, for dif- 
ferent values of the detuning Z?, with Uaa = 4, Uam = 12, 
and Plot — 2 (Fig. 2]). We can see that going from L = 1 
to L = 20 (equivalent to turning on the hopping param- 
eters ta and t„i) just leads to small differences at small 
g. For large g the hopping can be neglected, and results 
for i = 20 converge to those for L = 1. Nevertheless it 
is crucial to keep working with the full lattice problem 
instead of the one-site problem, since this is required to 
access global quantities such as the superfluid density. It 
is also the only way to get results for nearly-continuous 
values of ptot- 

A completely different behavior occurs when consider- 
ing a non-commensurate density, for instance ptot = 4/5 
(Fig. O. We consider here the case D < 2Uaa for sim- 
plicity. For this density, atoms can be placed on the lat- 
tice without increasing the interaction energy. The same 
holds for the molecules if D < 0. But for small g and 
-D > it is energetically more favorable to have atoms 
only, because 2 atoms have kinetic energy 4 times more 
negative than 1 molecule {2{—ta) = 4(— t^)). As a result 
the molecular density is vanishing for g — and D > 
and grows when turning on g, until reaching the optimal 
density for large g, pm = Ptot/^ = 1/5, in contrast to Fig. 



5 



1.5 



a- 1 




— p. 


D=6 


L=l 




D=6 


L=20 


Pn 


D=6 


L=l 


T T P„ 


D=6 


L=20 



P. 


D= 


10 


L=l 


O O P, 


D= 


10 


L=20 


P„i 


D= 


10 


L=l 


O O P„i 


D= 


10 


L=20 




FIG. 4: (Color online) The atomic and molecular densities as 
functions of the conversion parameter g and different values 
of the detuning D, for Uaa = 4, Uam = 12, and total density 
of particles ptot = 2. 



31 where for D < 6 the density pm decreases with increas- 
ing g. The atomic density follows the inverse behavior, 
starts for pa = ptot and converges to the optimal value. 
Pa = Ptot/^ = 2/5. 
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FIG. 5: (Color online) The atomic and molecular densities 
as functions of the conversion parameter g at different values 
of the detuning D. Uaa ~ 4, Uam = 12, and total density of 
particles ptot ~ 4/5. 

Having analyzed the system for two specific values of 
the total density, it is now interesting to perfom a scan of 
all values of ptot- Figure [6] shows the atomic and molecu- 
lar densities as functions of the total density ptot- At low 
filling, the particles are dilute and the on-site repulsion 
between atoms prevent double occupancies, so no binding 
between atoms can occur and the number of molecules 
remains zero for all values of D considered. Thus the 
atomic density increases linearly with the total density. 
As the filling increases, double occupancies occur leading 
to the creation of molecules, and decreasing the atomic 
density. Increasing the filling further leads to an "in- 



version of population" where the number of molecules is 
greater than the number of atoms. This inversion of pop- 
ulation is optimal at ptot — 2 for the chosen parameters 
because double atomic occupancies have an energy cost 
of 2Uaa — 8, whereas the creation of a molecule has an 
energy cost of D. Adding more particles to the system 
produces a saturation of molecules, and extra atoms just 
see a constant potential. 
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FIG. 6: (Color online) The densities of atoms and molecules 
as functions of the total density ptot at different values of the 
detuning D. Uaa = 4, Uam = 12, and g = 0.5. 

In order to identify incompressible phases, it is useful 
to look at behavior of ptot{p) for different values of D 
and g (Fig. [7|). Let us recall that the slope of this curve, 
dp/dp, is proportional to the isothermal compressibil- 
ity kt- Thus each horizontal plateau indicates an in- 
compressible Mott phase. This does not imply that this 
phase is insulating, as will be shown below. For D = — 1 
or D = and small conversion g = 0.5 one can iden- 
tify two incompressible phases by the presence of Mott 
plateaus at ptot = 2 and ptot = 3. For those parameters 
the usual Mott plateau occuring in pure bosonic systems 
at Ptot = 1 is absent. This is because extra particles can 
be added beyond ptot = 1 without the need of creating 
double occupancies, by converting atoms into molecules. 
For Ptot = 2, the phase is incompressible because any site 
is occupied by a molecule. Thus adding an extra atom 
requires the formation of an atom/molecule pair, which 
has an energy cost of Uam- For ptot = 3, each site is 
occupied with an atom/molecule pair, and adding extra 
atoms leads to double occupancies with energy costs of 
Uaa- Thus the phase is also incompressible. For D ~ 6 
and g = 0.5, we recover a Mott plateau at ptot — 1 be- 
cause creating a molecule has an energy cost of D that 
cannot be overcome by the associated negative kinetic 
and conversion energies. For large g however, the Mott 
plateaus at ptot — 1 and ptot — 3 disappear. Indeed, in 
this regime the conversions between atoms and molecules 
occur and overcome the energy cost of having two atoms 
on a single site, as well as the energy cost of creating a 
molecule. Thus extra atoms at ptot = 1 and ptot = 3 
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can go either into a molecule or doubly occupied sites, 
without changing the energy by a value greater than the 
finite-size-lattice gap, which vanishes in the thermody- 
namic limit. However ptot = 2 is still incompressible 
because any site is occupied either by two atoms or by a 
molecule. Thus an extra atom can go only on a site oc- 
cupied by a molecule, leading to an energy cost of Uam- 
Moreover a conversion process can no longer take place 
on this site, and the system has to pay the price of hav- 
ing a molecule all the time with the associated chemical 
potential D. This explains the large width of the cor- 
responding Mott plateau: approximately Z? + g + Uam- 
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FIG. 7; (Color online) The total density as a function of 
the chemical potential, and different values of D and g, for 
Uaa = 4, and Uam = f2. The slope of these curves is pro- 
portional to the isothermal compressibility, and horizontal 
plateaus indicate phases that are incompressible but not nec- 
essarily insulating (see text). 



We now study the potential superfluidity of the mix- 
ture by analysing the fluctuations of the atomic and 
molecular pseudo windings {W*'^) and {W^), and 
the correlated and anti-correlated windings (W^cor) ^"^^ 
{W^nt) (Fig- mil deflned in section IV. To discuss the 
results, it is useful to consider the corresponding curve 
in Fig. [Z](D = 6, 5 = 0.5; green curve). For ptot = 1 and 
Ptot = 3 all windings and pseudo windings vanish, show- 
ing that the system is frozen for those densities. The 
corresponding phases are Mott insulators. However for 
Ptot = 2 only the correlated winding vanishes, meaning 
that there is no global flow of particles, regardless of be- 
ing atoms or molecules. But individual species are flow- 
ing, each in the opposite direction of the other, leading to 
a large value of the anti-correlated winding. The phase is 
incompressible like a Mott insulator, but a supercurrent 
occurs for each of the species. We will refer to this phase 
as " Super-Mott" [l^ . We show in the following that this 
phase extends deep into the large-;? region of the phase 
diagram. 



FIG. 8; (Color online) The winding as a function of the filling 
for Uaa ~ 4, Uam = 12, Q = 0.5, aud D = 6. Error bars of the 
order of the symbol sizes. 



C. Phase diagrams 

Finally, by reproducing Fig[7| for different sets of pa- 
rameters g and D we are able to draw two phase dia- 
grams, one in the (/i, g)-j)\a,ne (Fig. [5]) and one in the 
(/i, I?)-plane (Fig. [10]). We can identify the three in- 
compressible phases discussed above, namely two Mott 
phases for ptot = 1 and ptot — 3, and the Super-Mott 
phase for ptot = 2. Those phases extend over regions of 
the phase diagram separated by superfluid regions. For 
small g, all incompressible phases are present. As g in- 
creases, the Super-Mott phase takes over the two Mott 
phases (Fig. For small or negative D the Super-Mott 
phase takes over the ptot = 1 Mott phase, whereas for 
large D it is the ptot = 3 Mott phase which yields to the 
Super-Mott phase (Fig. [TU)l . 




FIG. 9: (Color online) The phase diagram in the {p, g) plane, 
for Uaa = 4, Uam = 12, and D = 6. 
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FIG. 10: (Color online) The phase diagram in the (/i, -D) 
plane, for Uaa = 4, f/am = 12, and y = 0.5. 



the number of particles of the first species can lead to 
an inversion of population, resulting in the number of 
molecules greater than the number of atoms. In addi- 
tion to the usual superfluid and Mott phases occuring in 
boson models, we have identified an exotic "Super-Mott" 
phase, characterized by a vanishing compressibility and a 
superflow of both species but with anticorrelations such 
that there is no global supercurrent. Finally, we have 
produced two phase diagrams as a potential guide to 
detect the exotic Super-Mott phase. Since the Super- 
Mott phase occupies a big part of the phase diagrams, 
we expect it to be observable in experiments. We are 
currently investigating the model using a newly devel- 
oped algorithm that provides access to Green func- 
tions and momentum distribution functions, which can 
be measured in experiments. This will allow a direct 
comparison between theory and experiments. 



VI. SUMMARY AND DISCUSSION 

We have studied a two-species Bose-Hubbard model 
including a conversion term between the two species. 
Our model can be of interest for ultracold atom experi- 
ments using Feshbach resonances. The competition be- 
tween the kinetic, potential, and conversion terms leads 
to rich phase diagrams. We have shown that increasing 
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